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ABSTRACT 


The Double Pulsar, PSR JO737—3039A/B, has offered a wealth of gravitational experiments in the strong-field regime, all of which 
general relativity has passed with flying colours. In particular, among current gravity experiments that test photon propagation, the 
Double Pulsar probes the strongest spacetime curvature. Observations with MeerKAT and, in future, the Square Kilometre Array 
(SKA) can greatly improve the accuracy of current tests and facilitate tests of next-to-leading-order (NLO) contributions in both orbital 
motion and signal propagation.We present our timing analysis of new observations of PSR J0737—30394, made using the MeerKAT 
telescope over the last three years. The increased timing precision offered by MeerKAT yields a 2 times better measurement of Shapiro 
delay parameter s and improved mass measurements compared to previous studies. In addition, our results provide an independent 
confirmation of the NLO signal propagation effects and already surpass the previous measurement from 16-yr data by a factor of 1.65. 
These effects include the retardation effect due to the movement of the companion and the deflection of the signal by the gravitational 
field of the companion. We also investigate novel effects which are expected. For instance, we search for potential profile variations 
near superior conjunctions caused by shifts of the line-of-sight due to latitudinal signal deflection and find insignificant evidence with 
our current data. With simulations, we find that the latitudinal deflection delay is unlikely to be measured with timing because of 
its correlation with Shapiro delay. Furthermore, although it is currently not possible to detect the expected lensing correction to the 
Shapiro delay, our simulations suggest that this effect may be measured with the full SKA. Finally, we provide an improved analytical 
description for the signal propagation in the Double Pulsar system that meets the timing precision expected from future instruments 


such as the full SKA. 
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1. Introduction 


The Double Pulsar PSR JO737—3039A/B is a rich laboratory for 
strong-field gravity experiments. The system consists of a 23- 
ms recycled pulsar ("A") and a 2.8-s “normal” pulsar (“B”) in 
a nearly edge-on and slightly eccentric 2.45-hr orbit (Burgay 
et al. 2003; Lyne et al. 2004). Various relativistic effects have 
been precisely measured in previous works (Kramer et al. 2006b, 
20212), including periastron precession, time dilation (gravita- 
tional redshift and second-order Doppler effect), Shapiro delay 
due to light propagation in the curved spacetime of the compan- 
ion, and the orbital period decay which currently provides the 
most precise test of quadrupolar gravitational wave predicted by 


general relativity (GR). In addition, the relativistic spin preces- 
sion of B was measured by Breton et al. (2008) and the relativis- 
tic deformation of the orbit was newly detected in this system 
(Kramer et al. 2021a). All these make it a still unique system for 
gravity experiments. 

Comparing with other gravity experiments, the Shapiro de- 
lay measured in the Double Pulsar probes the strongest space- 
time curvature (~ 107?! cm?) in a precision experiment with 
photons, i.e. the interaction between gravitational and electro- 
magnetic fields (Wex & Kramer 2020). In addition, with 16 yr of 
data, Kramer et al. (20212) were able for the first time to measure 
higher-order effects of signal propagation in the strong gravita- 
tional field of a neutron star, which are currently not accessible 
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via any other method. These include retardation effect due to the 
movement of the companion (B) and aberrational light deflec- 
tion by the gravitation of the companion. The latter confirms the 
prograde rotation of A, which is consistent with the results mea- 
sured by Pol et al. (2018) using the emission properties of B and 
what is expected from binary evolution models. 

In this work, we present observations of PSR J0737—3039A 
with the MeerKAT telescope, a precursor for the Square Kilome- 
tre Array (SKA) at mid-frequency range. Thanks to its location 
in the Southern Hemisphere, it permits a timing precision more 
than two times better than that of the Green Bank Telescope for 
this pulsar (Bailes et al. 2020; Kramer et al. 2021b). This su- 
perior precision enables an independent and improved measure- 
ment of signal propagation effects within a very short time span. 
We also investigate effects that have been expected but not been 
studied in detail before. These include potential profile varia- 
tions due to latitudinal deflection, the detectability of latitudi- 
nal deflection delay, and the prospects of measuring the effect of 
lensing on the propagation time separately. 

This paper is organised as follows. Section 2 describes the 
MeerKAT observations on the Double Pulsar and data process- 
ing. In Section 3, we introduce the concepts of gravitational sig- 
nal propagation effects, including higher-order contributions to 
the Shapiro and aberration delay. The timing results and mass 
measurements are presented in Section 4. We then provide an 
in-depth study on the higher-order signal propagation effects in 
Section 5, with focus on latitudinal deflection and lensing. In 
addition, we provide an improved analytical description for the 
signal propagation in the Double Pulsar. Finally, we discuss the 
results and future prospects in Section 6. 


2. Observations and data processing 
2.1. MeerKAT observations 


The observations presented in this paper come from the 
MeerKAT telescope as part of the MeerTIME project (Bailes 
et al. 2020), which performs timing of known pulsars with var- 

ious scientific themes. Observations on PSR JO737—3039A are 
! conducted under the Relativistic Binary theme (RelBin, Kramer 
et al. 2021b), which focuses on testing the relativistic effects in 
binary pulsars to achieve measurements of neutron star masses 
and tests of theories of gravity. MeerTime observations are gen- 
erally recorded using the Pulsar Timing User Supplied Equip- 
ment (PTUSE) signal processor. This processor receives chan- 
nelised tied-array beamformed voltages from the correlator- 
beamformer engine of the MeerKAT observing system and is 
capable of producing coherently de-dispersed full-Stokes data in 
both filterbank (search) mode and fold (timing) mode, where the 
data are folded at the topocentric period of the pulsar. Details on 
pulsar observing set up with MeerTime are explained by Bailes 
et al. (2020). 

PSR J0737—30394 is regularly observed with a typical ca- 
dence of one month and duration of 3 hr. As the orbital period 
of this pulsar is ~ 2.45 hr, the observations are scheduled to start 
shortly before an eclipse and finish after the second eclipse, in 
order to observe the eclipses twice in one observing session. The 
session is typically composed of a 30-min observation with fold 
mode and search mode in parallel, followed by a 2-hr fold-mode 
observation and another 30-min fold-search dual-mode obser- 
vation. This specific arrangement is designed to maximise our 
sensitivity in detecting signal propagation effects, as well as in 
studying the magnetosphere of pulsar B (Lower et al. in prep.). 
Observations are performed with two receivers: The L-band re- 
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ceiver that covers the frequency range 856-1712 MHz and the 
UHF receiver that covers the frequency range 544-1088 MHZ, 
both with 1024 channels. The data presented here starts in March 
2019 and runs to May 2022. For the analysis in this paper, we 
use 29 full-orbit timing observations and 62 search-mode eclipse 
data sets, which amounts to a total of ~ 87 hr. 


2.2. Timing data reduction 


The raw 8 s-folded timing data from the PTUSE machines are 
processed with the MEERPIPE data reduction pipeline. MEERPIPE 
carries out radio frequency interference (RFI) removal using 
a modified coAsrGUARD algorithm (Lazarus et al. 2016), fol- 
lowed by flux and polarisation calibration. Details on polarisa- 
tion and flux calibration are described in Serylak et al. (2021) 
and Spiewak et al. (2022), respectively. 

After processing with MEERPIPE, the calibrated data prod- 
ucts are reduced using pulsar software package psRcHIVE! 
(Hotan et al. 2004). We first correct for the rotation mea- 
sure (RM) with the value measured in Kramer et al. (2021b), 
i.e. RM-120.84 rad m. As the L-band observations between 
March 2019 and February 2020 are restricted to 928 frequency 
channels (dropping 48 channels each from the top and bottom 
bands), to maintain consistency throughout the analysis, we re- 
duce the later L-band data to the same frequency channels. We 
treat the UHF-band data in the same way, as the roll-off ad- 
versely affects sensitivity of the top and bottom bands. 

For this system, a complete timing model is only available in 
the pulsar timing software Tempo? (Nice et al. 2015, more details 
will be given in Section 4). Therefore, to fold the data more ac- 
curately, all data are supplied with a polyco-format ephemeris 
with the values measured in Kramer et al. (2021a). Since the 
Double Pulsar rapidly changes its orbital phase, the time span 
(TSPAN) of a predicted pulse phase solution has to be as small 
as possible to keep a good precision?. With the PsRCHIVE ver- 
sion 2022-01-14, we set TSPAN to the minimum possible value 
which is 3 min. 

A known data processing issue with this pulsar is that the 
pulsar moves rapidly to a different orbital phase during the dis- 
persion delay time, hence the pulses received at the same time 
at different frequencies correspond to different orbital phases, 
therefore can not be folded with the same phase prediction. 
If not properly accounted for, this folding issue will cause 
frequency-dependent orbital smearing. Standard pulsar software 
like PsRcHIvE does not take this effect fully into consideration 
even with frequency-resolved rEMPO2 predictor^. To avoid this 
issue, we first de-disperse the total intensity data so that all fre- 
quencies correspond to the same orbital phase, then average the 
data first in frequency and then in time?. Because of the profile 
frequency evolution and scintillation effects, data are sub-banded 
in frequency, with 32 sub-bands for the UHF-band and 16 sub- 
bands for the L-band. 

As for the time-averaging, the integration time needs to be 
short enough in order to properly resolve the Shapiro delay and 


! http://psrchive.sourceforge.net/ 

2 http://tempo.sourceforge.net/ 

3 Our analysis suggests that the choice of TSPAN has a significant 
impact on the Shapiro parameters, a larger TSPAN leads to a large de- 
viation from the expected values. 

^ This issue is going to be addressed in psrcuive 2.0 under develop- 
ment. 

5 The order of processing matters. If reversed, the pulse phase appears 
to be different and phase offsets may be introduced. 
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Table 1. Information on MeerKAT observation and data set“ for PSR JO737—3039A. 


Receiver Centre Frequency Bandwidth Numberof Number of Time span Number of 
(MHz) (MHz) channels” sub-bands (MJD) TOAs 
L-band 1283.582 7715.75 928 16 58568 - 59721 83930 
UHF-band 815.734 493 928 32 58936 - 59663 137451 
“ Information presented here are for the trimmed data set, see Section 4.1. 
^ Effective usable channels. 
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Pulse phase 


Fig. 1. Pulse profile of PSR JO737—3039A observed at multiple fre- 
quencies with the MeerKAT UHF- and L-band receivers. 


the next-to-leading-order (NLO) signal propagation contribu- 
tions (quio. see Section 3), which are largest at superior con- 
junction. We perform a simulation to test the measurability of 
these NLO contributions with different integration time. The re- 
sults show that a good measurement of Shapiro delay and qur o 
can be achieved if the integration time is $ 32 s, but it becomes 
significantly worse if the integration time is longer than 1 min 
for quio and 2 min for Shapiro delay. Therefore, we average 
all data with 32 s integration time, consistent with the analysis 
by Kramer et al. (2021a). After frequency and time averaging, 
data are re-dispersed to allow measurement of dispersion mea- 
sure (DM) in the timing analysis. 


2.3. Wide-band templates and TOA extraction 


Wide-band observations like MeerTIME can suffer significant 
profile evolution in frequency, hence the traditional 1D tem- 
plate is not favoured. To best determine the pulse time-of- 
arrivals (TOAs) at multiple frequencies, we employ frequency- 
dependent 2D templates. With this technique, DM measure- 
ments are dependent on the DM value used to align (de-disperse) 
the 2D template. Due to the correlation between DM and pro- 
file evolution, DM measurements are to some extend frequency 
dependent, which can lead to a DM offset between L-band and 


is missing in our case. Therefore, to avoid this problem, we 
choose a bright observation from each band for making 2D tem- 
plates, and measure DM using data from their overlapping fre- 
quencies. Then, we use these DM values to de-disperse the cor- 
responding full-bandwidth data. This minimises the DM offset 
between L-band and UHF-band data which can be seen in Fig. 2. 
These data are then sub-banded and averaged in time. Finally, 
by smoothing the profiles with psrsmooth/PsRcHIVE, we obtain 
2D templates, with 16 sub-bands at L-band and 32 sub-bands at 
UHF-band. These templates are then used to measure frequency- 
resolved TOAs by cross-correlating with the reduced data using 
pat/psrcuive. The pulse profile of PSR JO737—3039A at mul- 
tiple frequencies is shown in Fig. 1. More information on the 
observing systems and data sets is given in Table 1. 


2.4. DM variation 


The wide-band observation and high precision of MeerKAT tele- 
scope make it possible to obtain an accurate DM measurement 
on a per-epoch basis so as to minimise the influence of DM noise 
in the data. To do so, we fit for only DM and spin frequency v 
for each observing epoch using 4-min TOAs, and keep the other 
parameters fixed. The DM measurements are shown in Fig. 2. 
Following Kramer et al. (2021a), we use a modified version of 
TEMPO for our timing analysis, which corrects dispersive delays 
for each TOA based on the exact DM measurement of that epoch. 

One should note that our data set does not show the appar- 
ent DM variation as a function of orbital phase as was seen in 
Ransom et al. (2004). It had been demonstrated that this effect 
occurs due to an unaccounted Doppler shift of the observational 
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frequency as the pulsar moves in a binary system®, this will be 
revisited by Hu, Porayko et al. (in prep.). More thorough inves- 
tigation of this effect as well as the frequency-dependent orbital 
smearing (see Section 2.2) is ongoing and will be presented in 
detail in the future publication. 


3. Signal propagation effects at superior 
conjunction 


In this section, we recapitulate the necessary concepts of signal 
propagation effects in the Double Pulsar, including the NLO con- 
tributions in the Shapiro delay and aberration delay, which were 
described in greater detail in Kramer et al. (2021a). 

Being a nearly edge-on binary system (i.e. i ~ 90°), the 
curved spacetime of the companion star (pulsar B) has a signifi- 
cant effect on the propagation of the pulsar’s signal. To leading- 
order this is the well-known Shapiro delay (Shapiro 1964), 
which is expressed in the following form for binary pulsars 
(Blandford & Teukolsky 1976; Damour & Deruelle 1986): 


ATO = —2rinAy, (1) 
A, =1 — er cos u — s [sin w (cosu — ez) 


+(1—24)'” cos wsinu] . (2) 


Here, u denotes the eccentric anomaly (from Kepler’s equation 
with eccentricity er), and w denotes the longitude of periastron 
measured from the ascending node. The time eccentricity er cor- 
responds to the eccentricity parameter in the Damour-Deruelle 
(DD) timing model (Damour & Deruelle 1986) that can be fitted 
in pulsar timing software TEMPO or TEMPO2 (Hobbs et al. 2006). 
The two post-Keplerian (PK) parameters r and s represent the 
range and shape of Shapiro delay, respectively. The shape pa- 
rameter is generally identified with the sine of the orbital incli- 
nation i as s = sini, whereas the range parameter is linked to 
the mass of the companion mp, which in GR follows r = Tg mg. 
The constant T = (GMN /c?, where c is the speed of light in 
vacuum and (GM) = 1.327 1244 x 107 cm? s7? is the nominal 
; solar mass parameter defined by the IAU 2015 Resolution B3 
(Prša et al. 2016). Through out the paper, all masses expressed 
in solar mass Mo are referred to the nominal solar mass by tak- 
ing the ratio Gm;/ (GMN (i = A, B), where G is the gravitational 
constant. 

The leading-order expression Eqs. (1) and (2) were obtained 
by integrating along a straight line (in harmonic coordinates) 
and assuming a static mass distribution when the pulsar’s signal 
propagates through the system (Blandford & Teukolsky 1976). 
In reality, the pulsar’s signal propagates along a curved path due 
to the deflection in the gravitational field of the companion and 
leads to a lensing correction to the Shapiro delay. This actually 
results in a reduced propagation time as a consequence of Fer- 
mat’s principle (Perlick 2004). The effect of lensing is not yet 
observable in any pulsar systems, but for completeness, one can 
extend Eq. (2) by an adapted version of the approximation in 
Klioner & Zschocke (2010, Eq. 73): A, — Ay + gA with 


5AM = 2rc[ag , (3) 


where the semi-major axis of the relative orbit ag = (x + xg)/s, 
with x and xg’ being the projected semi-major axes of pulsar A 
and pulsar B, respectively. For the Double Pulsar, one needs to 


é https://arxiv. org/e-print/astro-ph/0406321v2 
7 xp had been observed in Kramer et al. (2006b). 
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ascending node 


Fig. 3. Simplified illustration of effects related to the deflection of A’s 
radio signals (solid red) in the gravitational field of B (top down and side 
perspective). The observer is located at large distance along the x-axis. 
Apart from modifications in the propagation time due to a curved path 
in the gravitational field of B (lensing), one has a longitudinal deflec- 
tion delay (one) due to the fact that the pulsar has to rotate by more 
than 360° between two pulses while approaching the conjunction. Af- 
ter conjunction, it is less than 360°, which makes pulsar signals arrive 
earlier at the observer. In addition, there is a latitudinal effect, due to a 
latitudinal shift in the emission direction towards the observer. This can 
lead to changes in the pulse profile since the line of sight cuts a different 
part of the emission region, which can also be accompanied by changes 
in the pulse arrival times (more details in Sections 5.1 and 5.2). 


account for the fact that the companion star moves while the pul- 
sar’s signal propagates across the system. This effect is known 
as retardation effect or 1.5PN correction to the Shapiro delay 
(Kopeikin & Schäfer 1999; Rafikov & Lai 2006a). To sufficient 
approximation, the signal propagation delay can be extended to 


As = -2rIn(Ay + OAR? + SA3”), (4) 


where the retardation correction 5A‘ can be taken directly from 
Kopeikin & Scháfer (1999, Eq. 130) as 


2n x m . 2ms X mA 
aA coe Pe o es 
s Py mg (1 — eq)! 2 Py mg 
: 251/2 : 
sin w (cos u — er) + (1— eq) ^coscsinu 


(cosu — er) cos w — (1 — e2) 2 sin wsinu 
eq COS w + 
l—ercosu 


(5) 


The quantity Py denotes the orbital period, and ma denotes the 
mass of pulsar A. Note, in the Double Pulsar, the mass ratio 
ma/mg can be obtained in a theory-independent way (Kramer 
et al. 2006b; Damour 2007). Hence, apart from the Shapiro shape 
parameter s, Eq. (5) contains only Keplerian parameters. 
Moreover, the classical aberration expression (Smarr & 
Blandford 1976) assumes a flat spacetime for the propagation 
of the pulsar signals, which is no longer sufficient for describ- 
ing the observations of the Double Pulsar, particularly near the 
superior conjunction of pulsar A. One needs to account for the 
gravitational deflection of the pulsar’s signal caused by its com- 
panion (Doroshenko & Kopeikin 1995; Rafikov & Lai 2006b), 
which adds a lensing correction to the classical aberration. For 
pulsar A the misalignment angle between its spin vector and the 
orbital angular momentum is very small (< 3.2°, Ferdman et al. 
2008, 2013), which is in line with a low-kick birth event (cf. Pi- 
ran & Shaviv 2004; Willems & Kalogera 2004; Willems et al. 
2006; Stairs et al. 2006; Tauris et al. 2017). Since the spin of 
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A is practically parallel to the orbital angular momentum, the 
aberration delay can be simplified as 


AA — (sin V t ersin w) + p (6) 


The first term on the right-hand side of Eq. (6) is the classical 
aberration delay, where y = w + 0 is the longitude of pulsar 
with respect to the ascending node (0 is the true anomaly, which 
defines the angle between the direction of the pulsar and the pe- 
riastron), and the aberration coefficient 

x 


A= ~ 3.65 us. 7 
vPy(1 — e2)!/2 sin? i ü 29 


As A is practically not observable and can be absorbed by a 
shift in various timing parameters (see discussions in Damour 
& Deruelle 1986; Damour & Taylor 1992), we a priori add the 
aberration coefficient A as a fixed parameter in our timing model 
with the value given in Eq. (7). 

The second term pane in Eq. (6) is the higher-order correc- 
tion originating from the longitudinal deflection delay, and can 


be written as (Doroshenko & Kopeikin 1995) 


cos (yr + dy) l r 


londef _ : Em 
Ox =D A, + oARt , With p= 


TV X+Xp (8) 
Like in the Shapiro delay, retardation correction is also ac- 
counted for here. As a sufficiently good approximation, the po- 
sition of B when the signal reaches its minimum distance from 
B can be used (retardation corrected position; cf. Kopeikin & 
Schäfer 1999; Rafikov & Lai 2006a). The angle ó/*'* denotes 
the retardation related correction for the angle between the (co- 
ordinate) vector from B to A and the ascending node. 

As already discussed in Kramer et al. (20212), the NLO con- 
tributions in the Shapiro and aberration delays can not be tested 
separately in the Double Pulsar due to the similarity of their ef- 
fects on signal propagation. In addition, the lensing correction 
to the propagation delay (Eq. 3) is challenging to measure as it 
can be absorbed in the fit of Shapiro shape s (see Section 5.3 and 
discussions in Kramer et al. 20212). Therefore, to test the signifi- 
cance of the NLO contributions and to obtain an unbiased timing 
result, a common factor qwro is multiplied by these contributions 
and can be fitted for in our timing model: 


Art = Art x QNLO , 9) 
Ay" = AY" X dNLO » (10) 
goat a X QNLO - a 1) 


In GR, the scaling factor quio = 1. Figure 3 illustrates the dif- 
ferent effects related to signal deflection in the Double Pulsar 
system. 


4. Timing results 


For the timing analysis, we use the timing model in Tempo known 
as DDS, which is a modification of the DD model (Damour 
& Deruelle 1986) that uses a different parameterisation of the 
Shapiro delay. In DDS, the Shapiro shape parameter s is replaced 
by the logarithmic Shapiro shape parameter zs via 


Zs * —1n(l- s), (12) 


which is more suitable when s is very close to 1 (see Kramer 
et al. 2006a, 2021a). The NLO contributions in the Shapiro and 
aberration delays are also implemented in the latest DDS model, 
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Fig. 4. Post-fit residuals of PSR J0737—30394A using the DDS binary 
model as a function of time (top panel) and orbital phase of pulsar A 
with respect to the ascending node y (bottom panel). The MeerKAT L- 
band data are plotted in blue, whereas the UHF-band data are in red. 
The epochs of the excluded 6 observations are marked as black crosses. 


which can be measured through a common factor qw; o. Because 
the analytic inversion of the timing model developed in Damour 
& Deruelle (1986) is no longer sufficient for the Double Pulsar, 
primarily due to the NLO contributions, a numerical inversion of 
the timing model was also implemented in the latest DDS model 
in TEMPO (see Kramer et al. 20212). 


4.1. Timing parameters 


In our analysis, the full MeerKAT data set shows a large devi- 
ation in the Shapiro range parameter r compared to the 16-yr 
result (Kramer et al. 2021a). We perform a drop-out analysis by 
removing each observing epoch and fitting the parameters. We 
find that r is dependent on specific observing epochs, where 6 
epochs affect r by a significant amount while the rest epochs 
don't. These 6 epochs are marked as black crosses in Fig. 4. Af- 
ter excluding all these 6 epochs, r is consistent with the 16-yr 
result and the mass measurement in GR (Kramer et al. 20212). 
Even though a number of tests and simulations have been made, 
we are still unclear about the cause of this problem. Possible 
reasons could be systematic errors in the observations or fold- 
ing techniques. Note that all data were folded with TEMPO2 phase 
predictor during observations, which has shown outliers in this 
pulsar and has been doubly confirmed by our simulations. These 
outliers disappear after reinstalling a TEMPO polyco ephemeris in 
data processing (see Section 2.2), but we cannot rule out under- 
lying problems due to folding technique. The results shown here 
are based on data processed with polyco scheme?. Anyway, this 
issue should not affect the measurement of NLO signal propaga- 
tion effects, which is the main focus of this paper. Therefore, we 
leave this question to future studies. In the following analysis, 
we use a trimmed data set which excludes these 6 epochs. 


8 With the same set of observations, data processed with TEMPO polyco 

and rEMPO2 predictor show a noticeable difference (~ 3c) in the 
Shapiro parameters, where the result with polyco is closer to the 16-yr 
results and shows a smaller y?. 
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Table 2 and Fig. 4 present the results obtained from fitting the 
TEMPO DDS model to the trimmed MeerKAT data set. We fix the 
proper motion (Ue, us) and parallax 7x to the more precise val- 
ues determined from the 16-yr timing and VLBI measurement 
(see Kramer et al. 2021a). As the time span of our data is not 
sufficient to obtain a reliable measurement of the orbital period 
derivative Py and the orientation of the orbit (wo « 0?) is not 
at a favourable position for a precise measurement of the Ein- 
stein delay amplitude yg, we choose to fix these parameters to 
the more precise measurements from 16-yr data (Kramer et al. 
20212). Fixing the above parameters has no impact on the mea- 
surements of signal propagation effects and masses. The two PK 
parameters that describe the relativistic deformation of the or- 
bit, 6, and 6g (Damour & Deruelle 1986), are also held fixed 
at the GR value in our analysis, as ó, cannot be measured (see 
Kramer et al. 20212) and 6, is not yet measurable with the cur- 
rent MeerKAT data. 


The values shown in Table 2 are the result of 1000 Monte- 
Carlo (MC) runs, where in each run, a random realisation of 
proper motion, parallax, DM, yg, and È, is selected. The DM 
value is selected according to the DM measurements and un- 
certainties shown in Fig. 2. We use the aforementioned modified 
version of TEMPO to correct DM for each TOA and fit for all other 
timing parameters in each run. The numbers shown in Table 2 are 
the mean values of the distribution of each parameter after 1000 
MC runs, whereas the uncertainties is taken from the larger one 
among the standard deviation of the distributions and the maxi- 
mum error from TEMPO in all MC runs. 


In order to allow direct comparisons with previous publi- 
cations, parameters shown in Table 2 are measured with re- 
spect to the same epochs and the same terrestrial time standard 
UTC(NIST)? within the timescale “Barycentric Dynamical Time 
(TDB)" as implemented in Tempo. Even though TDB runs at 
a slower rate than the “Barycentric Coordinate Time (TCB)”, 
which was recommended by IAU 2006 Recolution B3!° (see 
also Soffel et al. 2003), this choice does not have any impact 
on the results presented in this paper (see discussions in Kramer 
et al. 2021a). To transfer the TOAs from topocentric to the Solar 
| System Barycentre (SSB), the DE436 Solar System ephemeris 
published by the Jet Propulsion Laboratory! is used. 


All binary parameters in Table 2 are consistent with the 16- 
yr data, except for x being different by ~ 30. This is because x 
is highly correlated with yg which is kept fixed in our fit. This 
should be improved in the future once we have enough MeerK AT 
data to fit for x and yg simultaneously. In our fit, the root mean 
square (RMS) is very close to the mean TOA uncertainty, and 
the reduced x? of individual observation is close to 1, suggesting 
that our result is not affected by jitter noise. We also perform 
simulations and single-pulse analysis following the methods in 
Parthasarathy et al. (2021) and find little evidence of jitter noise. 


The RMS of the MeerKAT data shown in Table 4 is more 
than two times better than that of the Green Bank Telescope (see 
Kramer et al. 2021b). Thanks to this much improved precision, 
the measurements of the Shapiro parameters improve quickly. 
Compared to Kramer et al. (2021a), the shape parameter s im- 
proves by a factor of 2 and the range parameter r improves by a 
factor of 1.3 (see Table 2). 


9 https://www.nist.gov/pml/time-and- frequency-division/ 
time-realization/utcnist-time-scale-Q. 

10 https://www.iau.org/static/resolutions/IAU2006_ 
Resol3.pdf 

!! https://ssd. jpl.nasa.gov/planets/eph export.html 
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Table 2. Timing parameters for PSR J0737—30394 using Tempo DDS 
binary model. Numbers in parentheses are 1o uncertainties referred to 
the last digits, obtained from the standard deviation of 1000 MC runs 
or maximum error from TEMPO, whichever is larger. The overall reduced 
X? is 0.99. 


Parameter Value 
Solar System ephemeris DE436 
Terrestrial time standard UTC(NIST) 
Timescale TDB 
Position epoch (MJD) 55045.0 
Timing epoch, to 55700.0 


Astrometric parameters 

Right ascension (R.A.), œ (J2000) 
Declination (Dec.), ó (J2000) 
Proper motion R.A., fg (mas yr^!) 
Proper motion Dec., us (mas yr^!) 
Parallax, 7x (mas) 


Spin parameters 

Rotational frequency (freq.), v (Hz) 
First freq. derivative, v (Hz s^!) 
Second freq. derivative, v (Hz s?) 


Binary parameters 

Orbital period, Py (days) 
Projected semi-major axis, x (s) 
Eccentricity, er 

Epoch of periastron, To (MJD) 
Longitude of periastron, w (deg) 
Periastron advance, & (deg yr!) 
Orbital period derivative (1072), P, 
Einstein delay amplitude, yg (ms) 
Logarithmic Shapiro shape, Zs 
Range of Shapiro delay, r (us) 
NLO factor for signal prop., gnLo 


Derived parameters 
s=sini=l-e* 

Orbital inclination, i (deg) 
Mass of pulsar A, ma(Mo) 
Mass of pulsar B, mp (Mo) 
Total mass, M (Mo) 


07:37:51.248 121(26) 
—30:39:40.705 36(42) 
—2.567(30)* 
2.082(38)* 
1.36(+0.12,—0.10)* 


44.054 068 642 001(56) 
—3.415 92(37) x 107 
—9.5(12) x 10727 


0.102251 559297 2(29) 
1.415 028 299(88) 
0.087 777 036(48) 

55700.233017 54(10) 
204.753 72(36) 
16.899 321(37) 
—1.247 920(78)* 
0.384 045(94)* 

9.669(77) 
6.163(16) 
0.999(79) 


0.999 936 9(+46/ -51) 
89.36(3) or 90.64(3) 
1.338 186(10) 
1.248 866(7) 
2.587 052(11) 


* Values adopted from Kramer et al. (20212). 


4.2. Mass measurements 


The standard approach for measuring the masses of a binary pul- 
sar system is using two PK parameters. Assuming GR, one can 
calculate the two a priori unknown masses based on the mea- 
surements of Keplerian parameters. For the Double Pulsar, the 
two most precisely measured PK parameters are periastron ad- 
vance w and the Shapiro shape parameter s. 

For the advance of periastron, in addition to the 1PN con- 
tribution, we also need to account for higher-order corrections 
due to 2PN effects and Lense-Thirring (LT) precession caused 
by spin-orbit coupling of pulsar A, as they are much larger than 
the measurement error of w (see Hu et al. 2020; Kramer et al. 
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Table 3. Mass measurements with a new modified DDGR model which 
accounts for NLO contributions in the orbital motion and signal prop- 
agation in this system. The MOI has been chosen to be J, = 1.28 x 
10?g cm? in the fit. 


Value 
1.338 186(10) 
1.248 886(5) 
2.587 050(8) 


Parameter 

Mass of pulsar A, ma (Mo) 
Mass of pulsar B, mg (Mo) 
Total mass, M (Mo) 


20212). For the analysis of this paper, the total intrinsic contri- 
bution to the periastron advance can be expressed, with sufficient 
precision, as (Damour & Scháfer 1988) 


PN :LLA 


w= OPN 4 PN qq, (13) 


The first and second post-Newtonian (PN) terms aN and w?PN 
are functions of masses and observed Keplerian parameters. 
While the situation is more complicated for the LT contribu- 
tion wT, as it requires the knowledge of the moment of inertia 
(MOD) of pulsar A, which is still not very constrained because of 
our limited knowledge of the equation of state (EOS) of neutron 
stars. As discussed in Hu et al. (2020) and Kramer et al. (2021a), 
one could measure the masses and the MOI simultaneously by 
introducing a third PK parameter Py into the @—s test. Such a test 
have already been made using the 16-yr data with an upper limit 
obtained: I, < 3.0 x 10? gcm? with 90% confidence (Kramer 
et al. 2021a). This measurement is expected to improve with the 
combination of the 16-yr data with MeerKAT data in a forthcom- 
ing paper, and should improve considerably over the next years 
with more data taken. This promises an important complemen- 
tary constraint on the EOS (Hu et al. 2020). For the calculations 
here, we take the value from Kramer et al. (2021a, Eq. 35) which 
uses the constraints on the EOS from Dietrich et al. (2020): 


c^ = —4.83(29, 35) x 107* deg yr! . (14) 


The Shapiro shape parameter s is the sine of the orbital in- 
, clination i. In Newtonian gravity, the orbital inclination is linked 
to the projected semi-major axis x via the binary mass function 
(e.g. Lorimer & Kramer 2004): 


(45) 


where x and the orbital frequency ny = 27z/Py are both mea- 
sured Keplerian parameters. Eq. (15) gets modified by a 1PN 
term in the 1PN approximation for Kepler's third law (see Eq. 3.7 
in Damour & Deruelle 1985 and Eq. 3.9 in Damour & Taylor 
1992): 


MPP 2/3 
sini = S | " (3 = mat) (ToMn) | 
TY mB 3M 


(16) 


Taking the measurements of Py, x, and masses from Table 2, 
one can calculate that the IPN correction is approximately 
1.27 x 1075. This correction was considered for the first time in 
pulsar analysis by Kramer et al. (20212), where the significance 
was about 1.30. Now with MeerKAT data, this 1PN correction 
is 2.50 significant and can not be ignored in the analysis. We 
hereby use the full IPN mass function Eq. (16) to measure the 
masses. 

Combining the PK parameters œ and s, one obtains the 
(Doppler-shifted) masses, which are listed in Table 2. These 


Residuals (us) 


L L L L L mn 
70 75 80 85 90 95 100 105 110 


y= w 0 (deg) 


Fig. 5. Aggregated residuals (blue) due to NLO contributions in Shapiro 
and aberration delay, shown in the orbital phase y. Residuals are re- 
scaled by (1 + er cos 0)! to account for secular variations in amplitude 
due to the precession of periastron. The black curve indicates the fitted 
qNLo (see Table 2) with the 20 range shown in grey shadow, which 
agrees very well with the theoretical prediction indicated by the red 
dotted line. 


measurements are fully consistent with those obtained with 16- 
yr data (Kramer et al. 20212), and the precision of ma and mg 
are also better. 

Alternatively, one can fit for masses using the timing model 
known as DDGR (Taylor & Weisberg 1989), which is based on 
the DD model where the PK parameters are explicitly calculated 
from the masses and the Keplerian parameters assuming GR. Be- 
side the Keplerian parameters, it fits explicitly for the total mass 
M and the companion mass mg. To make the measurements, 
we modify the DDGR model so that it incorporates all NLO 
contributions that need to be accounted for in this system, in- 
cluding NLO signal propagation effects, LT contribution & T^, 
NLO gravitational wave damping and mass loss contribution to 
Py (see Hu et al. 2020; Kramer et al. 2021a). An MOI needs to 
be provided to the model for the calculation of «^^ and the 
mass loss contribution to Py. For periastron advance c, the un- 
certainty from the MOI is still smaller than that from MeerK AT 
observations (see Eq. 14 and Table 2). Therefore, based on the 
EOS constraint from Dietrich et al. (2020), we fix the MOI to 
Ix = 1.28 x 10?g cm? in our fit. Table 3 shows the mass mea- 
surements obtained using the DDGR model. The results are fully 
consistent with the measurements derived from the DDS model, 
with smaller uncertainties in mp and M. 

Following Kramer et al. (20212), one could test the agree- 
ment of r in GR by comparing mt? = r/Ts = 1.251233) Mo 
(cf. Table 2) with the companion mass determined here, which 
gives 


7s [,88 = 1.0019(26). (17) 


This leads to a 5.3 x 107? (95% confidence) test of GR. 


5. Studying NLO signal propagation effects 


Because the Double Pulsar system is nearly edge-on to our line- 
of-sight (LOS, see i in Table 2), it is ideal for measuring signal 
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Fig. 6. Example of one eclipse observation at UHF band, plotted in 
intensity against orbital phase y and pulse phase. The intensity modula- 
tion occurs when pulsar A is eclipsed by the magnetosphere of pulsar B. 
Each integration is a sum of 8 pulses. When plotting, the discontinuities 
between recordings are patched with the previous sub-integration. 


propagation effects caused by the gravitational field of the com- 
panion near superior conjunction. The leading-order expression 
Eq. (1) is no longer sufficient to describe the signal propagation 
in the Double Pulsar. Such a model would result in significant 
residuals near superior conjunction when aggregating residuals 
in the orbital phase, as shown in Fig. 5. These residuals agree 
very well with the expected NLO contributions discussed in Sec- 
tion 3, which is shown by the red curve. The significance of the 
NLO corrections can be tested by scaling these corrections col- 
lectively by a common factor quio (cf. Eqs. 9-11; quy o-1 in 
GR) and fit for it. We find, with the trimmed data set, 


qno = 0.999(79) , (18) 
- which has surpassed the 16-yr result by 1.65 times with only 
~3 yr of data thanks to the much improved precision offered by 
MeerKAT. 

Following the definition of qw; o in Section 3, a fit for this pa- 
rameter involves two aspects of gravity: 1.5PN correction of the 
Shapiro delay due to the movement of the companion 6A‘*", and 
corrections related to the signal deflection in the gravitational 
field of the companion 6A!" and oo Even though these con- 
tributions can not be tested individually in a simultaneous fit be- 
cause of their similarity, one can still test one at a time while 
keeping the other one fixed (cf. Kramer et al. 2021a). We find 


qurLo[deflection] = 1.00(15) , 
qwrolretardation] = 1.00(17) . 


(19) 
(20) 


5.1. Searching for profile variation at eclipse 


The lensing correction to the aberration delay may not only lead 
to a shift in time in the longitudinal aspect but can also result 
in a change of the colatitude of the emission direction towards 
Earth, i.e., the latitudinal deflection delay. This would cause pro- 
file variations as the LOS cuts a different region of the pulsar 
beam (Rafikov & Lai 2006a,b). An illustration of the latitudinal 
deflection effect is shown on the right of Fig. 3. Various analyses 
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Fig. 7. Same as Fig. 6 but masking out the regions without pulses 
(blocked by the magnetosphere of pulse B). The red dashed lines in- 
dicate the orbital phase bins used in Figs. 8 and 9. 


have confirmed that pulsar A is an orthogonal rotator (Guille- 
mot et al. 2013; Ferdman et al. 2013; Kramer et al. 2021b), that 
the main pulse and the interpulse come from opposite magnetic 
poles. Therefore, we do not expect shifts of pulse components in 
phase as discussed in Rafikov & Lai (2006b) based on the (in- 
correct) assumption of an aligned rotator suggested by Demorest 
et al. (2004). 

The profile variation is expected to be maximum at the su- 
perior conjunction and symmetric around y — 90? (retardation 
corrected). This study requires high time resolution, for which 
we use the search mode data. We select the data that are near the 
eclipses and fold them into single pulses using TEMPO polyco 
(with TSPAN- 1 min). Data are then combined, cleaned, and po- 
larisation calibrated before integrated into total intensity and av- 
eraged in frequency. As the single pulses are still very weak, we 
average over every 8 pulses to increase the single-to-noise (S/N). 
An example of eclipse data is shown in Fig. 6. 

In order to get a high S/N profile, we first mask the regions 
where pulsar A's emission is blocked by pulsar B (see Fig. 7), 
and split the data in orbital phase with a step of Ay = 0.25? for 
89? < y < 91°. Then for each phase interval, we integrate pulses 
from all observations of a given band (L or UHF) together to 
increase the S/N. The resulting profiles are shown in Fig. 8 for 
the L-band data and in Fig. 9 for the UHF-band data, which are 
summed from 25 and 37 eclipses respectively. The difference 
between the added profiles at the eclipse and a 2-hr integrated 
profile at the non-eclipse part of the orbit is insignificant. The 
subtle residual structures in these figures can result from inter- 
stellar medium effects (DM variation and scintillation) based on 
our simulations. Therefore, we conclude that the current data are 
not (yet) sensitive to profile variations caused by the latitudinal 
aberration delay, or are not significant in the region that is seen 
by our LOS. These profiles will be used in a subsequent study 
on the geometry of the system. 


5.2. Simulation for latitudinal deflection delay 


To investigate whether the deflection delay caused by latitudinal 
deflection is measurable from pulsar timing, we perform a sim- 
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Fig. 8. Left: The black lines indicate the integrated profile of L-band data over various orbital phase intervals summed from 25 eclipses. The red 
curves indicate the reference profile integrated over a 2-hr observation excluding the eclipse part. The baseline of each profile is placed at the 
central orbital phase of the interval, and the number on the right side of the profile (np) indicates the upper estimate of the number of pulses in that 
interval. Right: Residuals of the added profile compared to the reference profile. 
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Fig. 9. Same as Fig. 8 but for the UHF-band data summed from 37 eclipses. 
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Fig. 10. Lensing simulation: pre-fit residual plotted against orbital phase 
w. Data displayed here are centred on w = 180° and span a decade. 
The scattering at the lower end of the curve is due to the precession of 
periastron. 


ulation based on a simple emission model, which consists of a 
| set of circular cones. Following Doroshenko & Kopeikin (1995) 
and Rafikov & Lai (2006b), the latitudinal deflection delay for 
pulsar A can be written as 


cos i sin (y + dy") 
(A, + ONT) tan yo í 


i oo eo (21) 


| where yo is the angle between the arc connecting the LOS 
and spin axis and the arc connecting LOS and magnetic axis. 
Note, Eq. (21) is based on the approximation of Doroshenko & 
| Kopeikin (1995) for the deflection angle, and therefore assumes 
, that the impact parameter is (sufficiently) large compared to the 
| Einstein radius. While this is sufficient, at least until full SKA 
' becomes operational, we present an improved description fur- 
ther below in Section 5.4. 

We include this deflection time delay in our test model as- 
suming yo = 45° !?, i.e. a relatively large latitudinal deflection 
delay, and scale it with a factor q!*““*f. We simulate high preci- 
sion TOAs using this test model and fit for the scaling factor. The 
pre-fit residuals show an advance signature with an amplitude of 
, —2.8 us and symmetric to y = 90°, which is in a similar shape 
to Shapiro delay but with opposite sign and smaller amplitude. 
However, after fitting for Shapiro parameters, the signature gets 
mostly absorbed and leaves residuals below 42 ns at superior 
conjunction. 


5.8. Prospects of lensing measurement 


Even though the retardation and deflection effects can be tested 
separately while keeping the other one fixed as shown in 
Eqs. (19) and (20), measuring the lensing correction to Shapiro 
delay SA!" independently is challenging. As already pointed out 
by Kramer et al. (20212), this effect is difficult to observe be- 
cause of its strong covariance with s, or equivalently z,. Our 
simulation also confirms that the lensing signature can be mostly 


12 Note that yo is not a constant, but our purpose here is to get a feeling 
for the measurability of the effect in timing rather than a proper account 
for the effect, which requires knowledge of the latitudinal variation in 
the emission pattern of the pulsar. Furthermore, as discussed at the be- 
ginning of Section 5.1, the beam geometry adopted by Rafikov & Lai 
(2006b) is not the correct one anyway. 
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Fig. 12. Uncertainty of factor q^" as a function of time span for the 
simulated data assumed to be 10 times better than the SKA 1. 


absorbed by z, in timing due to its symmetry with respect to con- 
junction. For demonstration purposes, we simulate 1 ns TOAs 
which model all NLO signal propagation contributions, then 
keep the retardation and deflection delay fixed in the model and 
fit for the scaling factor attached to the lensing correction, q'*" 
(corresponding to qwro in Eq. 10). Fig. 10 shows the residuals 
when lensing correction is not taken into account, leading to a 
reduced propagation time of about 850 ns as a result of Fermat's 
principle. Whereas after fitting for zs, the lensing signature gets 
absorbed and leaves the residuals to be below 70 ns (Fig. 11), 
making a detection with the current timing precision certainly 
impossible. 


To investigate whether lensing can be measured separately in 
the near future, we simulate TOAs for MeerK AT, MeerKAT ex- 
tension, and the first phase of the SKA (SKA 1) until 2030 based 
on similar assumptions made in Hu et al. (2020). In addition, 
as the TOA precision reduces significantly due to the intermit- 
tent signals during the eclipse (see Fig. 6), we account for this 
in our simulations by increasing the uncertainty of these TOAs 
based on MeerKAT observations. As a simple estimate, we as- 
sume GR and perform the simulation using the modified DDGR 
model with a grid fit to q'*^. If lensing is measurable, the value 
of q'*" should be close to 1. However, it turns out that with the 
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toward observer 


B’ 


Fig. 13. Schematic picture of the lensing geometry as used in Sec- 
tion 5.4. B’ denotes the retardation corrected position of B (cf. Klioner 
& Kopeikin (1992); Kopeikin & Schäfer (1999)). In principle there is 
a second photon path towards the observer (below B’). However, for 
the Double Pulsar this signal is not only significantly weaker, the path 
also comes so close to pulsar B that the photons are absorbed by the 
plasma-filled magnetosphere of B (cf. Lai & Rafikov 2005; Rafikov & 
Lai 2006b). 


observed and simulated data from 2019 to 2030, the uncertainty 
of q'*? is still larger than 1. 

To further push the precision, we assume that in the future 
an instrument will be available providing a timing accuracy one 
order of magnitude better than that of the SKA 1 (1.e., 100 ns for 
an integration time of 30 s), e.g., a future full SKA. As a rough 
estimate, here we only consider radiometer noise and ignore any 
other noise sources, such as jitter noise or scintillation noise. The 
uncertainty of q'" against the time span of the simulated data is 
shown in Fig. 12. With such precision, one would expect to get 
' a 5-c test of lensing with ~ 4 yr of data. 

From the simulation, we also obtain an estimated uncertainty 
for the common factor of NLO contributions qwro in the near 
future. Assuming no jitter noise, with MeerKAT and the SKA 1, 
we can expect a 50-o detection by 2030. 


5.4. Improvements in the timing model for < 50 ns precision 


Equations (8) and (21) are based on the approximation for the 
signal deflection used by Doroshenko & Kopeikin (1995). As 
discussed in detail in Kramer et al. (20212), this is still suffi- 
: cient for the analysis of current timing data. For that reason, the 
analysis in this paper is still based on Doroshenko & Kopeikin 
(1995), which (including corrections for retardation) is already 
part of the TEMPO distribution. However, in the future we can ex- 
pect to obtain a timing precision of better than ~ 50 ns near con- 
junction (+1°), so that an improved treatment of the deflection 
is required. In a series of papers, Rafikov and Lai have used the 
standard lensing equation to treat the signal propagation in the 
Double Pulsar near conjunction (Lai & Rafikov 2005; Rafikov 
& Lai 2006a,b). This allowed them to drop the assumption that 
the impact parameter is much larger than the Einstein radius. 
The standard lensing equation, however, is based on the assump- 
tion of small angles (see e.g. Schneider et al. 1992). Therefore, 
strictly speaking, Lai and Rafikov’s calculations are only valid 
near conjunction. Similar to calculations of Shapiro (1967) and 
Ward (1970), Wucknitz (2008) studied the deflection of photons 
in the gravitational field of a “point mass” for general lensing 
scenarios not limited to regions close to the optical axis. Based 
on these results, we give an analytic expression for the signal de- 
flection that is valid for the whole orbit, recovers the calculations 
by Lai and Rafikov near conjunction, and those of Doroshenko 
& Kopeikin (1995) if the impact parameter is large compared to 
the Einstein radius. 

In the following © is the angle between rap, and the direc- 
tion towards the observer, where rap denotes the vector from 


position of pulsar A at emission to the (retardation corrected) po- 
sition of pulsar B (the underlying geometry for our calculations 
is illustrated in Fig. 13). The deflection AO of A’s radio signal by 
pulsar B can be obtained from Eq. (24) in Wucknitz (2008), with 
the replacements a — AO, 0 — © + AO, and m > aż, where 
the quantity œg is the angle corresponding to the Einstein radius, 
and is given by 


2 {= 
p= - <l. 
€ V |rap] 


This is the maximum value AO can assume. The distance Dg 
in Wucknitz (2008) corresponds to Irag]. Consequently, one 
obtains 


(22) 


2 
q 
AO sin(O + AO) — - [1 + cos(O + Ag)] = O. (23) 
Since AO < ag « | forall angles ©, we can expand the equation 
above in AO while keeping terms only up to order a2. This leads 
to a quadratic equation 


2 


Q 

AG? + AO sin — EX 4 cos 9) « 0. (24) 

which has, under the assumptions made, one solution: 

A0 = 1 [ sin © + 202(1 + cos ©) - sino 25 
=> sin’ © + 2arc(1 + cos O) — sin : (25) 


For O « 1 this agrees with the standard lensing equation (see 
e.g. Schneider et al. 1992). 

The angle © € [0,7] needs to be computed from the 
retardation-corrected orbital phase via cosO = sinisin(y + 
ôy"). The longitudinal and latitudinal deflection delay are given 
by 


AO cos(y + dy") 
2zv  sinGsini 
AO sin + ov) 
~ Inv sin Otani tan yo ' 


respectively (cf. Eqs. (10) and (24) in Rafikov & Lai (2006b), 
with € = m — i and y = —7/2 spin of A aligned with orbital 
angular momentum). If © is much larger than ag one has AO ~ 
a2( 1 + cos ©)/(2 sin O). This corresponds to the approximation 
of Doroshenko & Kopeikin (1995) for the deflection angle, and 
one recovers Eqs. (8) and (21). 


t E (26) 


a = (27) 


6. Discussion 


In this paper, we have presented results from 3-yr timing ob- 
servations of the Double Pulsar using the MeerKAT telescope, 
with a specific focus on studying higher-order signal propaga- 
tion effects in the gravitational field of the companion. In order 
to minimise the effects from profile evolution and DM variation, 
we used frequency-dependent 2D templates to generate TOAs 
and a DM model to correct dispersive delay in TOAs. 

Thanks to its high inclination and orbital compactness, the 
Double Pulsar is a unique pulsar system for testing NLO signal 
propagation effects in strong fields. The significantly increased 
precision offered by MeerKAT permits an independent verifica- 
tion of NLO signal propagation effects and has already surpassed 


'3 Note, the Dg in Eq. (24) of Wucknitz (2008) is a typo and should not 
be there since it is already part of the definition of m. 
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the 16-yr result with only ~3 yr of data. In our analysis, the 
Shapiro shape parameter s has improved by 2 times compared 
to the previous result (Kramer et al. 2021a), which also leads 
to a better mass measurement. The Shapiro range parameter r 
agrees with GR at 5.3 x 107? (95% confidence). The precision 
of the measurement of NLO signal propagation effects qw; o has 
improved by 1.65 times. In this work, we investigated the poten- 
tial profile variation due to latitudinal deflection delay and the 
possibility of measuring lensing correction to the Shapiro delay, 
which has never been studied in detail before in pulsar analy- 
sis. With the current MeerKAT data, we found little evidence of 
profile variation at superior conjunction. It could be that the pro- 
file variation is not significant at the region we are looking at or 
our current data are not sensitive enough to identify it. We also 
did simulations on latitudinal deflection delay based on a simple 
emission model and found it unlikely to be detected because of 
its correlation with Shapiro delay. As for the lensing correction 
5A", we found it can be mostly absorbed by the Shapiro shape 
parameter. Our simulation showed that lensing is unlikely to be 
measured separately from timing before the full SKA or simi- 
larly powerful instruments, and may then be measurable with a 
few years of timing observations if noises like phase jitter and 
scintillation do not limit our precision. 

However, our analysis also showed that adding certain 
epochs has a significant impact on the measured Shapiro pa- 
rameters, but not on qwro. This could be due to the fact that 
the phase predicted using the polyco scheme is in particular 
worse at the superior conjunction, which caused the discrepan- 
cies in the Shapiro parameters. Comparison of polyco with dif- 
ferent TSPAN values showed residuals oscillating near the supe- 
rior conjunction, and we may already be limited by the precision 
of polyco scheme. Of course, there may exist other unknown 
systematic errors in the data. 

To support our timing analysis and study of profile variations 
at eclipse due to latitudinal deflection, we also checked profiles 
from all observations. We found variations in the total profile 
from epoch to epoch. The differences in the profiles are more 
prominent at lower frequencies and broadband. Our simulation 
suggested that these profile variations are likely to be associated 
with DM variation and scintillation. Even though we have sub- 
banded data into 16/32 frequency bands and used 2D templates, 
profile variations may still have an impact on timing. The study 
of profile variations will be continued in a subsequent work to 
improve the constraint on the geometry of the system. 

Moreover, although not discussed in the paper, we found ev- 
idence of red noise in the spectrum of timing residuals with an 
amplitude two orders of magnitude larger than for typical mil- 
lisecond pulsars. If not taken into account, it may strongly af- 
fect astrometric parameters, as well as influence binary parame- 
ters, according to our simulations. This makes it more difficult to 
combine the 3-yr MeerKAT data set and the 16-yr data set. Given 
that the timing precision of the former significantly outperforms 
the latter, the weighting of MeerKAT data already exceeds the 
16-yr data and can dominate noise modelling. For the purpose 
of this paper, we did not include 16-yr data because of their mi- 
nor contribution to the quy o measurement (~10% improvement). 
But for studying secular relativistic effects, an appropriate noise 
modelling may be required to combine these data. We investigate 
this in further ongoing studies. 

In the future, continuing observations with MeerKAT and the 
SKA will further improve the precision of tests on signal prop- 
agation effects, and a 50-o detection of qur o can be expected 
by 2030. For that, we have also provided an improved analytical 
description of the signal propagation in the Double Pulsar. Fur- 
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Table 4. Comparison of the MeerKAT timing precision Crys assumed 
in Hu et al. (2020) and from real observations with L-band and UHF- 
band receivers, scaled to 5-min integration time over the full bandwidth. 


Telescope / receiver Reference C'RMs (us) 
MeerKAT L band Hu et al. (2020) 1.06 
MeerKAT L band this work 0.90 
MeerKAT UHF band this work 0.55 


thermore, as demonstrated by Hu et al. (2020), the precision of 
secular relativistic effects will also be greatly improved and will 
eventually allow the measurement of the MOI of pulsar A and 
the NLO gravitational wave damping in the near future. The tim- 
ing precision of the MeerKAT data used in this work is even bet- 
ter than that assumed in Hu et al. (2020), which is based on early 
L-band data from MeerKAT (See Table 4). This makes their pre- 
dictions conservative and we are likely to achieve even better 
measurements with future observations. 
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